﻿#include "discrete_face_sphere.hxx"
#include "discrete_face_loop.hxx"
#include "idx_mesh_utils.hxx"
#include "acis/sphere.hxx"
struct scan_inter {
    double u;
    int low_inter;
    int up_inter;
};

void sphere_faceter::init() {
    // sphere is only closed on v , not on u
    // mainly [-pi/2,pi/2]*[-pi, pi)
    localSphere = (sphere*)&f->geometry()->equation_for_update();
    SPAinterval u_interval = f->geometry()->equation().param_range_u();
    SPAinterval v_interval = f->geometry()->equation().param_range_v();
    u_le = u_interval.start_pt();  // u的左边界
    u_ri = u_interval.end_pt();    // u的右边界
    v_le = v_interval.start_pt();  // v的左边界
    v_ri = v_interval.end_pt();    // v的右边界

    u_tolerance = 1e15;                                                                                // u上容差
    v_tolerance = 1e15;                                                                                // v上容差
    if (f->geometry()->equation().closed_u()) u_tolerance = localSphere->param_range_u().length() / 2;  // 判断跨边界所用值
    if (f->geometry()->equation().closed_v()) v_tolerance = localSphere->param_range_v().length() / 2;

    le_low.x = u_le;  // 左下角点
    le_low.y = v_le;
    ri_high.x = u_ri;  // 右上角点
    ri_high.y = v_ri;

    SPAunit_vector zeroNormal = f->geometry()->equation().eval_normal(SPApar_pos(0, 0));
    if (f->sense()) zeroNormal = -zeroNormal;
    if (zeroNormal != localSphere->uv_oridir) normalR = REVERSED;
}

void sphere_faceter::decideUVlen() {
    // 优化牛顿法，直接求取du和dv
    // 用参数域上(t,t),(-t,t),(t,-t)和(-t,-t)四个点来进行计算
    // 详细推导过程可以看对应的文档
    double r = localSphere->radius;  // 球面的半径

    // 计算角度容差
    float cos_nt = cos((nt) / 180 * M_PI);               // 将角度容差转为cos值，便于后续计算
    float t_nt_ab = asin(sqrt(2 - 2 * cos_nt)) / 2;      // 点A和点B满足角度容差的t
    float t_nt_ac = acos(sqrt((1 + cos_nt) / 2));        // 点A和点C满足角度容差的t
    float t_nt_ad = acos(sqrt(sqrt((1 + cos_nt) / 2)));  // 点A和点D满足角度容差的t
    // 计算面容差
    float t_st = acos(sqrt(1 - st / r));  // 满足面容差的t

    du = 2 * std::min({ t_nt_ab, t_nt_ac, t_nt_ad, t_st });  // 需要从上面四个t中选择一个最小的t来作为du，du=2t
    assert(!is_equal(du, 0.0));
    ulen = std::ceil((u_ri - u_le) / du);  // 因为要保证du和dv能整除Pi，所以需要计算ulen和vlen，再获得最终的du和dv

    vlen = ulen * 2;  // 这里要保证du和dv是一个正方形，所以vlen应该是ulen的两倍

    // 求得最终的du和dv
    assert((ulen != 0) && (vlen != 0));
    du = (u_ri - u_le) / ulen;
    dv = (v_ri - v_le) / vlen;
}

void sphere_faceter::attachMesh() {
    REVBIT r = f->sense();
    SPApar_pos param;

    // 面上采点
    for (double v = v_le; v < v_ri + SPAresabs / 2; v += dv) {
        // 构造一条从y = v 的直线 求交
        // 求交结果记录在inters中，以两个点为一组，两个点中间的点为inside的点，添加到points中
        SPAposition a(u_le - 2 * SPAresabs, v, 0.0);
        SPAunit_vector dir(1.0, 0.0, 0.0);
        straight v_line(a, dir);
        SPAinterval line_range(0.0, u_ri - u_le + 4 * SPAresabs);
        v_line.limit(line_range);
        std::vector<scan_inter> inter_u;

        for (unsigned i = 0; i < edges.size(); i++) {
            point curr = points[edges[i].v1()];
            point next = points[edges[i].v2()];
            if (curr == next) continue;
            SPAposition b(curr.x, curr.y, 0.0);
            SPAposition c(next.x, next.y, 0.0);

            SPAvector edge_dir = c - b;
            straight edge_line(b, normalise(edge_dir));
            SPAinterval edge_range(0.0, (edge_dir).len());
            edge_line.limit(edge_range);
            curve_curve_int* results_head = int_cur_cur(v_line, edge_line);
            curve_curve_int* results = results_head;

            if (results != nullptr) {
                if (results->next != nullptr) {
                    delete_curve_curve_ints(results_head);
                    continue;
                }
                double p1 = results->param1;
                SPAposition inters = v_line.eval_position(p1);
                unsigned j = 0;
                for (; j < inter_u.size(); j++) {
                    if (inter_u[j].u == inters.x()) {
                        break;
                    }
                }
                if (j == inter_u.size()) {
                    inter_u.push_back(scan_inter(inters.x(), 0, 0));
                }
                if (is_equal(v, curr.y) && is_equal(v, next.y)) {
                    inter_u[j].low_inter += 1;
                    inter_u[j].up_inter += 1;
                }
                else if (is_equal(v, curr.y) && v > next.y || v > curr.y && is_equal(v, next.y)) {
                    inter_u[j].low_inter += 1;
                }
                else if (is_equal(v, curr.y) && v < next.y || v < curr.y && is_equal(v, next.y)) {
                    inter_u[j].up_inter += 1;
                }
                else {
                    inter_u[j].low_inter += 1;
                    inter_u[j].up_inter += 1;
                }
            }

            delete_curve_curve_ints(results_head);
        }
        std::sort(inter_u.begin(), inter_u.end(), [](const scan_inter& a, const scan_inter& b) { return a.u < b.u; });
        for (unsigned i = 0; i < inter_u.size(); i++) {
            if ((inter_u[i].up_inter + inter_u[i].low_inter) % 2 == 1) {
                if (i == 0 || i == inter_u.size() - 1) {
                    inter_u.erase(inter_u.begin() + i);
                    i -= 1;
                    continue;
                }
            }
            if (inter_u[i].low_inter % 2 == 0 && inter_u[i].up_inter % 2 == 0) {
                inter_u.erase(inter_u.begin() + i);
                i -= 1;
            }
        }
        for (int i = 0; i < inter_u.size(); i += 2) {
            double inf = inter_u[i].u;
            if (i + 1 >= inter_u.size()) continue;
            double sup = inter_u[i + 1].u;
            for (double u = get_cloest_boundary_point(u_le, du, inf); u < sup; u += du) {
                if (u < inf) continue;
                points.push_back(point(u, v));
            }
        }
    }

    // 调用CDT
    CDT::Triangulation<double>* pcdt = nullptr;  // CDT进行三角剖分
    int _count = 0;

    if (pcdt) ACIS_DELETE pcdt;
    pcdt = ACIS_NEW CDT::Triangulation<double>;
    CDT::RemoveDuplicatesAndRemapEdges(points, edges);

    pcdt->insertVertices(points);
    pcdt->insertEdges(edges);
    pcdt->eraseOuterTrianglesAndHoles();

    CDT::TriangleVec triangles = pcdt->triangles;
    std::vector<CDT::V2d<double>> vertice = pcdt->vertices;
    INDEXED_MESH* mesh = ACIS_NEW INDEXED_MESH(vertice.size(), pcdt->triangles.size(), pcdt->triangles.size() * 3);
    for (int i = 0; i < vertice.size(); i++) {
        CDT::V2d<double> pp = vertice[i];
        SPApar_pos param(pp.x, pp.y);
        SPAposition pos = f->geometry()->equation().eval_position(param);
        SPAunit_vector normal = f->geometry()->equation().eval_normal(param);
        if (r) normal = -normal;
        mesh->add_vertex(pos, normal, param);
    }

    for (int i = 0; i < triangles.size(); i++) {
        int p1 = triangles[i].vertices[0];
        int p2 = triangles[i].vertices[1];
        int p3 = triangles[i].vertices[2];
        mesh->add_polygon(i, 3);
        indexed_polygon* poly0 = mesh->get_polygon(i);
        poly0->set_vertex(0, &mesh->get_vertex(p1));
        poly0->set_vertex(1, &mesh->get_vertex(p2));
        poly0->set_vertex(2, &mesh->get_vertex(p3));
    }
    af_update(static_cast<ENTITY*>(f), IDX_MESH_APP, IDX_MESH_ID, mesh);
}

outcome sphere_faceter::facet() {
    init();
    decideUVlen();

    auto generator = [&](SPApar_pos begin, SPApar_pos end) {
        double min_x = std::min(begin.u, end.u);
        double max_x = std::max(begin.u, end.u);
        double min_y = std::min(begin.v, end.v);
        double max_y = std::max(begin.v, end.v);
        bool is_u = begin.u == end.u;
        bool is_less = is_u ? (begin.v < end.v) : (begin.u < end.u);
        std::vector<SPApar_pos> tmp;
        if (is_u) {
            for (double point_y = begin.v; is_less ? point_y < end.v : point_y > end.v; point_y += is_less ? dv : -dv) tmp.emplace_back(begin.u, point_y);
        }
        else {
            for (double point_x = begin.u; is_less ? point_x < end.u : point_x > end.u; point_x += is_less ? du : -du) tmp.emplace_back(point_x, begin.v);
        }
        return tmp;
        };
    // Getloops处理
    auto loop_processor = LoopProcessor(f);
    loop_processor.RegisterGenerator(generator);
    loop_processor.Execute();
    for (auto x : loop_processor.points) {
        points.emplace_back(x.u, x.v);
    }
    for (int loop_idx = 0; loop_idx < loop_processor.loops_indexs.size(); loop_idx++) {
        auto& point_indexs = loop_processor.loops_indexs[loop_idx];
        int point_indexs_size = point_indexs.size();
        if (point_indexs_size < 2) continue;
        for (int i = 1; i < point_indexs_size; i++) {
            edges.emplace_back(point_indexs[i], point_indexs[i - 1]);
        }
        edges.emplace_back(point_indexs.front(), point_indexs.back());
    }
    attachMesh();
    return outcome();
}

outcome facet_face_sphere(FACE* f, double nt, double st) {
    sphere_faceter faceter(f, nt, st);
    return faceter.facet();
}
